body .main-container {
  max-width: 100% !important;
  width: 100% !important;
}
body {
  max-width: 100% !important;
}

body, td {
   font-size: 16px;
}
code.r{
  font-size: 14px;
}
pre {
  font-size: 14px
}

Let us set some global options for all code chunks in this document.

knitr::opts_chunk$set(
  message = FALSE,    # Disable messages printed by R code chunks
  warning = FALSE,    # Disable warnings printed by R code chunks
  echo = TRUE,        # Show R code within code chunks in output
  include = TRUE,     # Include both R code and its results in output
  eval = TRUE,       # Evaluate R code chunks
  cache = FALSE,       # Enable caching of R code chunks for faster rendering
  fig.align = "center",
  out.width = "100%",
  retina = 2,
  error = TRUE
)

rm(list = ls())
set.seed(1982)

1 Preprocessing

Let us now load some required libraries.

# Load required libraries

# inla.upgrade(testing = TRUE)
# remotes::install_github("inlabru-org/inlabru", ref = "devel")
# remotes::install_github("davidbolin/rspde", ref = "devel")
# remotes::install_github("davidbolin/metricgraph", ref = "devel")
library(INLA)
library(inlabru)
library(rSPDE)
library(MetricGraph)

library(plotly)
library(dplyr)
library(tidyr)
library(MASS)
library(glmnet)
library(car)
library(ggplot2)
library(plotly)

library(sf)

library(here)

We now define some R functions with self-explanatory names.

# Functions to manipulate the data
standardize <- function(vector) {return((vector - mean(vector)) / sd(vector))}
minmax_scaling <- function(x) {return((x - min(x)) / (max(x) - min(x)))}
convert_to_binary <- function(input_vector) {return(ifelse(input_vector != 0, 1, 0))}

1.0.1 Prepare the data

The following chunk uses the function creates_covariates_on_mesh() defined on Preprocessing/28Window.case.file4.addcovariatesonmesh.R to build a mesh (on the same graph as the one used in this file) and compute the value of all available covariates on the mesh nodes.

Go to Preprocessing 8.0.1 Raw code to see the body of the creates_covariates_on_mesh() function.

h = 0.05
# source(here("Preprocessing/28Window.case.file4.addcovariatesonmesh.R"))
# creates_covariates_on_mesh(h)
# Sys.sleep(5)

################################################################################
################################# PREPARE THE DATA #############################
################################################################################

# loading the data
load(here("Data_files/data_on_graph_with_covariates_no_consecutive_zeros_partialtomtom.RData"))
load(here("Graph_objects/graph_construction_28_03_2024partialtomtomwhichlonglatsf.RData"))
load(here("Data_files/data_on_mesh_with_covariates_partialtomtom.RData"))


data = data_on_graph_with_covariates %>%
  mutate(across(starts_with(c("class_", "upto")), list(ind = convert_to_binary))) %>%
  mutate(across(c("bus", "signal", "stop", "crossing"), ~round(., 5))) %>%
  mutate(across(c("SpeedLimit", "density_per_hour"), standardize)) %>% dplyr::select(-datetime)


mesh = data_on_mesh_with_covariates %>%
  mutate(across(starts_with(c("class_", "upto")), list(ind = convert_to_binary))) %>% # this creates new columns
  mutate(across(c("bus", "signal", "stop", "crossing"), ~round(., 5))) %>%
  mutate(across(c("SpeedLimit", "density_per_hour"), standardize))

1.0.2 Add data and build mesh

Let us add the speed observations to the graph and build the mesh. Observe that mesh already has the value of the available covariates on the mesh nodes.

sf_graph$add_observations(data = data, group = "day", normalized = TRUE, clear_obs = TRUE)
sf_graph$build_mesh(h = h)
Sys.sleep(5)

2 Modeling

2.0.1 Stationary model

  • Observe that we are considering replicates.
stat.time.ini <- Sys.time()
################################################################################
################################# STATIONARY MODEL #############################
################################################################################

rspde_model_stat <- rspde.metric_graph(sf_graph,
                                         parameterization = "matern",
                                         nu = 0.5)

data_rspde_bru_stat <- graph_data_rspde(rspde_model_stat,
                                        repl = ".all",
                                        name = "field")

st.dat.rep.stat <- inla.stack(
    data = data_rspde_bru_stat[["data"]]["speed"], 
    A = list(data_rspde_bru_stat[["basis"]],1), 
    tag = "est",
    effects =
      list(c(data_rspde_bru_stat[["index"]], list(Intercept = 1)), 
           list(SpeedLimit = data_rspde_bru_stat[["data"]]["SpeedLimit"])
      )
    )


f_stat = speed ~ -1 + Intercept + SpeedLimit + f(field, model = rspde_model_stat, replicate = field.repl)

rspde_fit_stat <-
  inla(f_stat,
      data = inla.stack.data(st.dat.rep.stat),
      family = "gaussian",
      control.predictor = list(A = inla.stack.A(st.dat.rep.stat))
  )

stat.time.fin <- Sys.time()
print(stat.time.fin - stat.time.ini)
## Time difference of 15.23668 secs
summary(rspde_fit_stat)
## 
## Call:
##    c("inla.core(formula = formula, family = family, contrasts = contrasts, 
##    ", " data = data, quantiles = quantiles, E = E, offset = offset, ", " 
##    scale = scale, weights = weights, Ntrials = Ntrials, strata = strata, 
##    ", " lp.scale = lp.scale, link.covariates = link.covariates, verbose = 
##    verbose, ", " lincomb = lincomb, selection = selection, control.compute 
##    = control.compute, ", " control.predictor = control.predictor, 
##    control.family = control.family, ", " control.inla = control.inla, 
##    control.fixed = control.fixed, ", " control.mode = control.mode, 
##    control.expert = control.expert, ", " control.hazard = control.hazard, 
##    control.lincomb = control.lincomb, ", " control.update = 
##    control.update, control.lp.scale = control.lp.scale, ", " 
##    control.pardiso = control.pardiso, only.hyperparam = only.hyperparam, 
##    ", " inla.call = inla.call, inla.arg = inla.arg, num.threads = 
##    num.threads, ", " keep = keep, working.directory = working.directory, 
##    silent = silent, ", " inla.mode = inla.mode, safe = FALSE, debug = 
##    debug, .parent.frame = .parent.frame)" ) 
## Time used:
##     Pre = 0.446, Running = 5.8, Post = 0.98, Total = 7.22 
## Fixed effects:
##              mean    sd 0.025quant 0.5quant 0.975quant   mode kld
## Intercept  20.883 0.374     20.163   20.877     21.638 20.877   0
## SpeedLimit  3.060 0.280      2.517    3.061      3.597  3.062   0
## 
## Random effects:
##   Name     Model
##     field CGeneric
## 
## Model hyperparameters:
##                                           mean    sd 0.025quant 0.5quant
## Precision for the Gaussian observations  0.012 0.000      0.012    0.012
## Theta1 for field                         2.853 0.070      2.717    2.852
## Theta2 for field                        -1.089 0.217     -1.508   -1.091
##                                         0.975quant   mode
## Precision for the Gaussian observations      0.012  0.012
## Theta1 for field                             2.993  2.849
## Theta2 for field                            -0.656 -1.102
## 
## Marginal log-Likelihood:  -55592.19 
##  is computed 
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')
fit.rspde = rspde.result(rspde_fit_stat, "field", rspde_model_stat)
summary(fit.rspde)

2.0.2 Nonstationary model

  • Observe that we are using the computed parameters from the stationary model as initial values for the nonstationary models.
nonstat.time.ini <- Sys.time()
################################################################################
############################# NON STATIONARY MODEL #############################
################################################################################

B.sigma = cbind(0, 1, 0, mesh$SpeedLimit, 0)
B.range = cbind(0, 0, 1, 0, mesh$SpeedLimit)

rspde_model_nonstat <- rspde.metric_graph(sf_graph,
                                            start.theta = c(fit.rspde$summary.log.std.dev$mode,
                                                            fit.rspde$summary.log.range$mode,
                                                            rep(0, 2)),
                                            B.sigma = B.sigma,
                                            B.range = B.range,
                                            parameterization = "matern",
                                            nu = 0.5)

data_rspde_bru_nonstat <- graph_data_rspde(rspde_model_nonstat,
                                        repl = ".all",
                                        name = "field")

st.dat.rep.nonstat <- inla.stack(
    data = data_rspde_bru_nonstat[["data"]]["speed"], 
    A = list(data_rspde_bru_nonstat[["basis"]],1), 
    tag = "est",
    effects =
      list(c(data_rspde_bru_nonstat[["index"]], list(Intercept = 1)), 
           list(SpeedLimit = data_rspde_bru_nonstat[["data"]]["SpeedLimit"])
      )
    )


f_nonstat = speed ~ -1 + Intercept + SpeedLimit + f(field, model = rspde_model_nonstat, replicate = field.repl)

rspde_fit_nonstat <-
  inla(f_nonstat,
      data = inla.stack.data(st.dat.rep.nonstat),
      family = "gaussian",
      control.predictor = list(A = inla.stack.A(st.dat.rep.nonstat))
  )

nonstat.time.fin <- Sys.time()
print(nonstat.time.fin - nonstat.time.ini)
## Time difference of 19.44462 secs
summary(rspde_fit_nonstat)
## 
## Call:
##    c("inla.core(formula = formula, family = family, contrasts = contrasts, 
##    ", " data = data, quantiles = quantiles, E = E, offset = offset, ", " 
##    scale = scale, weights = weights, Ntrials = Ntrials, strata = strata, 
##    ", " lp.scale = lp.scale, link.covariates = link.covariates, verbose = 
##    verbose, ", " lincomb = lincomb, selection = selection, control.compute 
##    = control.compute, ", " control.predictor = control.predictor, 
##    control.family = control.family, ", " control.inla = control.inla, 
##    control.fixed = control.fixed, ", " control.mode = control.mode, 
##    control.expert = control.expert, ", " control.hazard = control.hazard, 
##    control.lincomb = control.lincomb, ", " control.update = 
##    control.update, control.lp.scale = control.lp.scale, ", " 
##    control.pardiso = control.pardiso, only.hyperparam = only.hyperparam, 
##    ", " inla.call = inla.call, inla.arg = inla.arg, num.threads = 
##    num.threads, ", " keep = keep, working.directory = working.directory, 
##    silent = silent, ", " inla.mode = inla.mode, safe = FALSE, debug = 
##    debug, .parent.frame = .parent.frame)" ) 
## Time used:
##     Pre = 0.2, Running = 12.4, Post = 1.02, Total = 13.6 
## Fixed effects:
##              mean    sd 0.025quant 0.5quant 0.975quant   mode kld
## Intercept  20.799 0.405     20.005   20.796     21.604 20.796   0
## SpeedLimit  2.863 0.273      2.332    2.862      3.404  2.863   0
## 
## Random effects:
##   Name     Model
##     field CGeneric
## 
## Model hyperparameters:
##                                           mean    sd 0.025quant 0.5quant
## Precision for the Gaussian observations  0.012 0.000      0.012    0.012
## Theta1 for field                         3.011 0.114      2.784    3.012
## Theta2 for field                        -0.722 0.292     -1.302   -0.720
## Theta3 for field                         0.371 0.119      0.133    0.372
## Theta4 for field                         0.741 0.232      0.279    0.743
##                                         0.975quant   mode
## Precision for the Gaussian observations      0.012  0.012
## Theta1 for field                             3.234  3.015
## Theta2 for field                            -0.152 -0.713
## Theta3 for field                             0.602  0.375
## Theta4 for field                             1.193  0.751
## 
## Marginal log-Likelihood:  -55598.78 
##  is computed 
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')
summary(rspde.result(rspde_fit_nonstat, "field", rspde_model_nonstat))
save.image(here(paste0("Models_output/", rmarkdown::metadata$title, ".RData")))